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I. INTRODUCTION 



Nonlinear metamaterials (NMM) are rapidly developing fields extending the concepts of metamaterials as well as 
nonlinear optics. 

NMM are interesting because of their extraordinary properties compared with natural nonlinear materials, including, 
amongothers, enhanced nonlinearity due to resonance l|,l2|, novel phase-matching condition due to negative refractive 
index 3j |^ . Generally, NMM are inherently nonlinear lattices driven by external light radiation with subwavelength 
lattice constants. Under weak radiation, the responses of the lattice should be approximately described by the linear 
optics with appropriate modifications coming from nonlinearity. In this regime, low-order nonlinear coefficients in 
Taylor expansions as for natural nonlinear materials would be an appropriate description. Many recent papers work 
in this weak radiation regime 3|, [5|-|9|. When radiation becomes stronger but not strong enough to damage the 



building materials of meta-atoms, in some cases 2|. |10| there may exist a scale of radiation intensity in/above which 
the lattice has such complex responses due to modulational instability that it can hardly be described by uniform 
optical constants (e.g., linear and/or nonlinear refractive index). That scale is characterized by the possible existence 



of bistability under uniform assumption regardless of the modulational instability for the case in 2 
numerical simulations 
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Full-wave 



12] of the structures in a bistability proposal [13J and the results for a similar nonlinear system 
14] support that the bistability based on uniform assumption may be disturbed by modulational instability. In this 
paper we will study modulational instability for a different kind of NMM, magneto-elastic metamaterials (MEM). 

Recently a new kind of metamaterials called MEM are proposed [1 51 ] presenting nonlinearity arising from mechanic 
interactions among individual meta-atoms. In the sections below, we will hrst establish an effective circuit dynamic 
model for MEM and then investigate the properties, especially modulational instability, of four cases: (i) a dimer of 
split-ring resonators (SRR); (ii) a one-dimension (ID) dimcrizcd array or an array of dimers; (iii) a ID uniform array 
of SRR or a ID polymer; (iv) layers two-dimension (2D) SRR array. We find that except for the first case all the 
other three many-body cases can present modulational instability above an radiation intensity threshold, which can 
inevitably disturb the bistability (if any) under uniform assumption. 



II. THEORETICAL MODEL AND AN SRR DIMER 



A simple model based on effective circuit is developed in this section for an SRR dimer and generalized for the 
other many-body cases in the following sections. Linear stability analysis indicates that modulational instability can 
not occur for an SRR dimer. 



As proposed in 15[, an MEM is composed of SRR which can deviate from their original positions due to the 
magnetic dipole forces from their neighbors but subject to respective mechanic restoring forces. The consequence 
of the changing of positions is that the mutual inductance among neighbor SRR can thus be affected, rendering a 
feedback mechanism of the electromagnetic response of MEM, another way of saying nonlinearity. For the simplest 
case, two interacting SRR sharing the same central axis, we model the SRR as effective LCR circuits and further 
assume that the mechanic resonance due to the restoring force and the magnetic resonance of the LCR circuit are 
decoupled, i.e., the frequency of the latter is much larger than that of the former, and also that the damping of the 
mechanic oscillation is very fast compared to the varying rate of envelope of magnetic response. In other words, the 
balanced position of each SRR can simply be determined by the averaged magnetic forces (related to the envelope of 
magnetic response) it experiences and the restoring force coefficients. Possible electric interactions are not taken into 



account in this paper for simplicity. These assumptions are consistent with the ones made in 
here for deriving simple dynamic equations for the system. 

In Appendix we obtain the dynamic equations for a dimer under slowly- varying approximation 



15] and are intended 



2ifl 



dQi- 
dr 



{2.1 = Ul,2, 



(1) 



where Q, f2, r, 7, k and u are dimensionless quantities related to charge, frequency, time, damping, mutual inductance 
and electromotive force, respectively. Included in Appendix lAl are more details on the approximation to derived Eq.([T]) 
and the definitions of the dimensionless quantities. 
The stationary states of Eq. (|A8|) then satisfy 



{-n 2 + i-fn+l)Q h2 - kQ 2 1 + Q 2 Re(QiQ^) Q 2 ,i=«i,2 



(2) 



Under uniform radiation u± = 112 = u, the uniform stationary states Q\ = Q2 = Q satisfy 



(-Q 2 + ijtl + l)Q - nfl 2 l + n 2 \Q\ 2 Q = u 



(3) 



Bistability can easily be obtained from the above uniform equation. See Fig[l]for an example. 



II. 1. Modulational instability for an SRR dimer 

Will the two SRR always respond uniformly under uniform radiation u\ = U2? In order to answer this question 
we should investigate the modulational instability of the uniform stationary states in Eq.Q. Before we do linear 
stability analysis for the uniform equation ([3]) , we want to point out that the stationary equations ([2]) without making 
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FIG. 1. A typical bistability curve for the uniform response of a dimer. Parameters: Q — 0.73, 7 — 0.01, k = 0.8 



uniform assumption do not have nonuniform solutions 



Iff, implying no modulational instability could occur for a 



dimer. Linear stability analysis in Appendix iBl also confirms this speculation. 



III. ID DIMERIZED SRR ARRAY 

Since a dimer can only have uniform response, it will be convenient to treat it as a unit nonlinear element rather 
than two SRR. A question is: can modulational instability occur when an array of dimers are arranged along a line 
sharing the same central axis (or plane), i.e., for a ID dimerized SRR array which models the experimental structures 



in Fig. 5 of la] . In this section, we investigate the modulational instability of such an array. 

If each dimer in such an array is treated as a unit nonlinear element and the interactions between neighbors are 
considered to be linear, the dynamic equation can be easily written according to Eq.([T]), 

2i(l + re)fi^ + {-n 2 + ijQ + 1 - aft 2 - Kfl 4 \Q n \ 2 )Q n - K rf fi 2 (Q„_ 1 + Q B+1 ) = u. (4) 
dr 

where na is the linear interaction coefficient between nearest dimers. We have denoted Q by Q for simplicity and 
employed the nearest-neighbor approximation. 

Under uniform assumption and periodic boundary conditions, the stationary uniform states can be obtained from 

{-fl 2 + + kQ 2 - nn 4 \Q\ 2 )Q - 2n d Q 2 Q = u. (5) 



Bistability can also present for Eq.© as shown in Fig 12 
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FIG. 2. Typical uniform responses of a dimerized array. Parameters: 7 = 0.01, k = 0.8. The dashed portions of the curves are 
instable ranges. 

III.l. Modulational instability of a dimerized array 



We do stability analysis in Appendix [C] for this case. The growth rate for each Fourier component of the fluctuation 
is calculated. Unlike the dimer, a dimerized array can indeed present modulational instability which may disturb the 
upper or lower branch of the bistability depending on the sign of the interaction coefficient Kd (A minus Kd means 
that the dimers share the same central plane). The unstable ranges are denoted by the dashed line in Figj2] 

In order to see to what extent modulational instability can disturb the uniform response, we do time-domain simu- 
lations by directly integrating Eq.Q. In order to trigger the modulational instability, we introduce tiny fluctuations 
to the initial states of the array, i.e., Q,i(t = 0) randomly distribute within [0, 10~ 6 ]. The radiation is turned on 
at t — and is stationary thereafter. We consider an array containing nine dimers and monitor the time-domain 
responses of the third, sixth and ninth one in Fig|3]and FigUJ 

In FigOJ the simulations are for the Kd = 0.2, = 0.660 case in FigJ5] The upper branch of the bistability is 
completely covered by the unstable range. We test two radiation intensities corresponding to u = 0.006 (corresponding 
to a stable uniform state outside the unstable range) and u = 0.01 (inside the unstable range). For u = 0.01, the 
modulational instability manifests itself as chaos. 

In FigUl the simulations are for the Kd = —0.2, f2 = 0.80 case in FigJ5J Part of the lower branch of the bistability 
is unstable as shown in FigJ^l We test two radiation intensities corresponding to u = 0.01 and u = 0.02. For the 
latter, we test two initial conditions, Q n (r — 0) ~ and Q„(r = 0) « 0.22. For u = 0.01, no modulational instability 
exists as we expected from FigJ5J However, when u = 0.02 and Q n (r = 0) « 0, the system without experiencing the 
possible nonuniform responses directly skip to the upper branch which does not suffer from modulational instability. 
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FIG. 3. Time-domain responses of three dimers in the dimer array. Parameters: k — 0.8, 7 = 0.01, Hd = 0.2, Q = 0.660, 
u = 0.006 (a), u = 0.01 (b). 



If we prepare the initial states inside the unstable range Q n (r = 0) « 0.22 with tiny fluctuations, we can indeed 
observe nonuniform responses due to modulational instability for a while, but the system will finally leak into the 
upper uniform branch anyway. 

u=0.01 b) u=0.02, Q (t=0)=0 C ) 




u=0.02, Q (t=0)=0.22 



real(Q 6 ) 
-real(CL) 



0.5 



-0.5 



real(Q 6 ) 
-real(Q 3 ) 
-real(QJ 







real(Q 6 ) 
— real(Q 3 ) 
— real(Q 9 ) 







1000 2000 3000 4000 5000 



1000 2000 3000 4000 5000 



2000 4000 6000 8000 



FIG. 4. Time-domain responses of three dimers in the dimer array. Parameters: k = 0.8, 7 = 0.01, Hd = —0.2, fl — 0.80, 
u — 0.01 (a), u = 0.02 (b and c). The initial conditions are Q n (r = 0) ~ (a and b) and Q n (r = 0) m 0.22 (c) 



IV. ID UNIFORM ARRAY 

In the ID array in the previous section, the nonlinear interactions occur pairwise. In this section, we consider the 
case in which the nonlinear interaction is considered to be uniform, i.e., each SRR has equal nonlinear interactions 
with its neighbors. More precisely, the balanced position of each SRR is determined by the magnetic forces of all its 
neighbors (rather than only that of its partner as in the dimerized array case) and the positions of SRR in turn affect 
the mutual inductances. In Appendix Q3J we derive the dynamic equations for this case under the slowly- varying and 
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nearest-neighbor approximations, 

2iti^- + {-n 2 + in 7 + i)Q n = u n + K n 2 [i + n 2 Re(Q n Q* n+1 - 2Q„_xQ; + Qn-iQn_ 2 )]Qn-i 

+ K Q 2 [1 + n 2 Re{Q n+1 Q* n+2 - 2Q n Q* n+1 + Q n Q* n -i)]Qn+i- (6) 
The stationary states under uniform assumption and periodic boundary conditions are 

[-(1 + 2k)^ 2 + 1]Q = u. (7) 

We observe that the relation between the response Q and the radiation u is linear and thus no bistability can exist. 
That is because the array is infinitely long (due to the periodic boundary conditions) so that the net magnetic force 
vanishes for each SRR. But it is not the case for an array of finite length, e.g., a dimer in SecHH or for a nonuniform 
array, e.g., a dimerized array in Sec lIIIl 

While the uniform responses do not present nonlinear effects, we can still study possible modulational instability 
of the uniform linear responses. That is analytically done in Appendix [D] in which we have obtained the growth rate 
for each Fourier component of the fluctuation. Only when all the modes of fluctuations decay, the uniform response 
is stable. In FigO we find the boundary between the stable and unstable range for k = 0.5 and 7 = 0.1. We can see 
that there exists a critical frequency (about Q = 0.74) below which no modulational instability occurs. 
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FIG. 5. Boundary between stable and unstable range for ID uniform array. 

Next, we do time-domain simulations to show how modulational instability manifest itself in this case. In Fig|6l 
we show three cases, one in the stable range, the other two in the unstable range. We find that the system can have 
nonuniform but stationary states in the unstable range for those specific cases. 
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FIG. 6. Time-domain responses of three SRR of a uniform array containing ten SRR with periodic boundary conditions. 
Parameters: k = 0.4, 7 — 0.1. 

V. 2D LAYERED SRR ARRAY 



151 ] (Fig.l). We construct a similar model for this 2D 



We then consider the 2D structures originally proposed in 
layered structures in Appendix [E] as the cases in previous sections. The main simplifications are nearest-neighbor 
approximations, slowly-varying approximations, and Taylor approximations for mutual inductances and magnetic 
forces. The the main deviation of these approximations from the model in |15| is that the former gives simple 
bistability of uniform responses as shown below while the the latter gives a more complicated one. But our model can 
capture the discrete nature of that structure which is essential for studying modulational instability. The calculations 
below show that modulational instability may also occur in the uniform bistability range just as in the dimerized-array 
case in SeclIID 

Based on the simplifications mentioned above, we derive dynamic equations for this 2D nonlinear structure in 
Appendix lEl Bistable responses under uniform assumptions can occur in this model (see FigJT]). But when doing 
time-domain integrations in FigjSJ we find that modulational instability manifest itself as nonuniform (e.g., limit- 
cycle) responses disturbing the uniform bistability. We also find that despite the complex responses within one of the 
two layers, the responses of counterparts in the two layers keep close to each other in the tested case, similar to the 
dimer in Secllll 



VI. DISCUSSION 



We discuss why modulational instability seems to usually accompany the uniform bistability for nonlinear many- 
body systems in this paper. Uniform bistability is characterized by a middle unstable branch, which according to the 
Fourier analysis in Appendix [B] is the unstable range for the uniform (or k = 0) fluctuations. Since for a infinitely 
large system the wave vector of fluctuations continuously ranges from k = to k = it /a (a is the lattice constant), 
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FIG. 7. Bistability under uniform assumption for 2D layered magneto-elastic metamaterials. Parameters: g = 0.5, mo = —0.1, 
mio = 0.5, m2o=0.1, 77121 = 0.1, to = 0.78, 7 = 0.01. The detailed definitions of the dimensionless parameters can be found in 
Appendix [E] 

it is natural that the whole unstable range, the union of unstable ranges for all wave vectors, can be obtained by 
continuously extending the original unstable middle branch. It can also be seen from the Fourier analysis in Appendix 
E] that the range of bistability disturbing by modulational instability is proportional to the coupling coefficients 
between neighbor elements. Therefore, for negligible coupling, modulational instability should not be an important 
factor. 



It is worth pointing out that a recent publication [14,] which studies the modulational instability accompanying 
possible uniform bistability also verify the arguments above in a completely different system. 



We develop an effective circuit model for a recently designed magneto-elastic metamaterials and study the modu- 
lational instability of them. Four cases are studied: a dimer, ID dimerized array, ID uniform array, and 2D layered 
array. We find that a dimer do not present modulational instability and can be treated as a unit nonlinear element. 
All the other three many-body cases can have modulational instability showing chaotic, stationary but nonuniform, 
or limit-cycle responses. The modulational instability seems to usually accompany the uniform bistability for the ID 
uniform array and 2D layered array. 
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VII. CONCLUSION 
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FIG. 8. Time-domain responses of three SRR of a 2D layered array containing 2x10x10 SRR with periodic boundary conditions 
in the directions parallel to the two layers. The parameters are the same as in Fig[7] Two intensities with modulational 
instability are tested, and limit-cycle responses are observed. The figure (c) shows the difference between two counterparts at 
the site (2,9). 

Appendix A: Dynamic equations for an SRR dimer 



The dynamic equations for the two SRR based on effective circuit model are 



L 



dt 2 



R- 



dt 



C 



= £t 



Ni- 



di 2 



(Al) 



where q\^. is the respective charge hold by the capacitors of the two SRR, L is the self-inductance, R is the resistance, 
C is the capacitance, M is the mutual inductance, and £ e tuJt is the effective electromotive force. 

Since the mutual inductance is determined by the averaged magnetic force between SRR, we are only interested in 
the amplitudes Qi : 2 of the charges qi^ = Qi,2{t)e. lu)t . Then we should derive the dynamic equations for <3i,2 under 
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slowly- varying approximation as in 10], which finally read 



2iujL^hl + (-uj 2 L + iloR + i)Qi >a - Mw 2 Q 2 ,i = f i,a 



(A2) 



Next, we remember that the mutual inductance M is a function of the magnetic response. First, we approximate 
M by its Taylor expansion reserving only the first two orders, 



MnMo + ^(b-b ), 



(A3) 



where b is the distance between the two SRR, 60 is the distance for zero magnetic force, and Mq is the mutual 
inductance for b = 6 - 

According to Hooke's law, 



K{b-bo) = -<F(b,h,I 2 )> 



(A4) 



where K is the coefficient of the restoring force, F is the instantaneous magnetic force as a function of the two currents 
7i 2 i and < • • • > denotes the time-averaged value. 

The instantaneous magnetic force should be of the form, 



F{b,h,I 2 )=hI 2 f{b)*hI 2 f{b ), 



(A5) 



while the time-averaged value should thus be 



<F(b,h,I 2 ) >« -Re(/iiJ)/ 



(A6) 



where 1\ 2 are complex amplitudes of the currents and fo = f(bo). Note that we have made further simplification 
that f(b) is approximated by f(bo). 

Then after combining Eq. (|A3IIA"6)) the mutual inductance M can be expressed as a function of the complex currents 

dM /oRe(JiiJ) 



M = Mi 







db 2K 



(A7) 



At this point we are ready to write down the full dimensionless dynamic equations for the SRR dimer under the 
slowly-varying approximation and other simplifications mentioned above, 



2ifl- 



dr 



+ (-fT +i-ffl + 1)Q 1)2 - l + fi 2 Re(QiQ; 



12,1 = "1,2, 



(A8) 



where we have defined \I C \ 2 = 

"1,2 = C£l,2Uo/\Ic\- 



2KM 

r DM 
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Appendix B: Linear stability analysis for an SRR dimer 

We do linear expansion of the dynamic equations (|TJ) near the uniform stationary states obtained from Eq.@. 
Substituting Q±. 2 — Q + 6i t2 to Eq.([T]) where Q is a uniform stationary state, we arrive at the linear dynamic 
equations for fluctuations 5\ t 2 

2in— — + {-n 2 + j 7 o + i)(5i, 2 - K n 2 (i + n 2 \Q\ 2 )s 2 1 

UT 

- k 9_(5*q2 + S 1 \Q\ 2 + S 2 \Q\ 2 + SIQ 2 )=0. (Bl) 
Subtract one of the two equations above with the other, 

2in d<y5l ~ 52) + f-o 2 + ijn + i){s 1 - s 2 ) - nn 2 (i + n 2 \Q\ 2 )(s 2 - s t ) = o. (B2) 

UT 

Clearly, this equation only has decaying solution for Si — 8 2 , which means that the uniform states from Eq.(lO) are 
stable under nonuniform fluctuations. 



Appendix C: Linear stability analysis for a dimerized array 

From the dynamic equations the dynamics of fluctuations around uniform stationary states Q n (t) — Q s + 5 n (t) 
should read. 

2i (i + K )n^ + (d + 2C 2 \Q s \ 2 )6 n + C 2 Q 2 s 8* n + C7 3 (<5„-i + S n+1 ) = (CI) 
where d = -fl 2 + i-yQ + 1 - nQ 2 , C 2 = -nfl 4 , C 3 = -K d Vl 2 . Or 



2i(l + k)Q-^ + AS n + B5* + C(S n -i + 5 n+1 ) = 



(C2) 



The above equation is linear and has periodic symmetry, so the eigenmodes should be uniform plane waves. Let 
S n = ae lkn + Be~ lkn . Then 



— = -irj[(A + 2C cos k)a + BB*] 
dr 

dB* 

-— = -ir)[-B*a - (A* + 2C cos k)B*] 
dr 

where l/i-q = 2i(l + «)fi. Note that 77 < 0. 

The eigen values of the coefficient matrix in the right-hand side are 



(C3) 
(C4) 



A = 7, 

= V 



7 f7± y/\B\ 2 - |A + 2C*cosfc| 2 +7 2 ^ 2 

7O ± V(^ 4 |Q S | 2 ) 2 - (-^ 2 + 1 - nil 2 - 2nn 4 \Q s \ 2 - 2K d Q 2 cosfc) 5 



(C5) 
(C6) 
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The A > indicates unstable states , which lie in the range 



^\Q S \ 2 G [~ (2a - Va 2 - 36 2 ) , ~ (2a + ^ a 2 - 36 2 



(C7) 



where a = —SI 2 + 1 — k£1 2 — 2njVl 2 cos k, b — jCt. 



Appendix D: ID uniform array and modulational instability 



We consider a periodic chain of SRRs with nonlinear magneto-elastic interactions and periodic boundary conditions. 
The deviations of each SRR from the original balanced positions are denoted by Xi, i — 0, 1, • • •. According to 
Hooke's law, and the nearest-neighbor approximation, 



Kx n = -Re(I„J* +1 - i" n _iI*)/ . 



We have further approximated the magnetic force by its zeroth order. 
Mutual inductance between the n'th and n+l'th is 

M(x n ,x n+ i) = M + ^-{x n+1 - x v 
ox 

Then the stationary equation is 



-oj LQ n + iujRQ Tl 



Qn 

~c 



£ 



dx 2K 



u 2 Q n ~i 



LU 2 Q n 



+ 1- 



(Dl) 



(D2) 



dx 2K 

where we have assumed nearest-neighbor approximation for the mutual inductance. 

The above equation can be recast into a dimensionless form in a similar way to the dimer case in Appendix fXJ 



(D3) 



or 



(-S1 2 + ifi7 + l)Q n = u n + kC1 2 [1 + Re(i n i* n+1 ~ 2i n - X i* n + i n -ii* n _ 2 )}Qn-i 

+Ktt 2 [l + Re(i„+ii* +2 - 2i n i* n+1 + i n i* n _ 1 )]Qn+i (D4) 



(-n 2 + ifi 7 + 1)Q„ = u n + K n 2 [l + fl 2 Re(Q n Q* n+1 - 2Q n ^Q* n + Q„_iQ;_ 2 )]Q„-i 

+ K n 2 [l + n 2 Re(Q n+1 Q: i+2 - 2Q n Q* n+1 + Q n Q* n - X )]Qn+i (D5) 



since i n = iilQ. 
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The dynamic equations can be easily obtained from Eq. (|D5|) under slowly- varying approximation, 

2t n^- + (-ft 2 + zft 7 + i)Q n = u n + K n 2 [i + n 2 Rc(Q n Q* n+1 - 2Q„_ 1 q; + Q n -iQ»_ 2 )]Qn~i 

CLT 

+ K ft 2 [l + n 2 Re{Q n+1 Q* n+2 - 2Q n Q* n+1 + Q n Q* n -i)]Qn+i (D6) 
In order to study the modulational instability, dynamics for fluctuations, Q = Q s + S n is also needed. 

2ift— 2 = -(-ft 2 + zft 7 + l)8 n + Kfi 2 (<5„_i + 5„+i) + K ft 2 ([*] + [*])Q (D7) 

UT 

where 

[*] = ft 2 Re[Q* (<5„_ 2 - 5 n _x - <J n + S n+1 )} (D8) 

and 

[4] = ft 2 Re[Q* (5 n _ x - *„ - <5„ +1 + <J„ +2 )] (D9) 

Then 

2ift— - = -(-ft 2 + ?:ft 7 + l)S n + Kft 2 (6„_i + 5 n+1 ) + (D10) 

UT 

K ft 4 Re[Q*(<5„_ 2 - 2<5„ + S n+2 )}Q (Dll) 

Let (5„ = ae ife ™ + ^e" ifcn . We have 

2iftd = -(-ft 2 + i-fVL + 1)q + 2ftft 2 acosfc + ftft 4 [Q*(cos2fc - l)a + Q(cos2fc - 1)/3*]Q (D12) 
2iSl/3 = -(-ft 2 + i 7 ft + l)(3 + 2kQ 2 f3 cosk + K ft 4 [Q* (cos2fc - l)/3 + Q(cos2fc - l)a*]Q (D13) 

We take the complex conjugation of the second equation 

2iftd = -(-ft 2 + i 7 ft + l)a + 2Kft 2 acosfc + K ft 4 [Q*(cos2fc - l)a + Q(cos2fc - 1)0*]Q (D14) 

-2ift/3* = -(-ft 2 -i 7 ft + l)/3* + 2/d7 2 /r cos fc + reft 4 [Q (cos 2/fc- l)/3* + Q*(cos2fc - l)a]Q* (D15) 

Reform 

a = — {Aa + B/3*) (D16) 

4= — (-B*a-A*/3*) (D17) 

where i?/ = 2ift (rf > 0), A = -(-ft 2 + i 7 ft + 1) + 2 K ft 2 cosfc + K ft 4 |Q| 2 (cos 2fe - 1) and B = K ft 4 Q 2 (cos2fc - 1). 
The characteristic polynomial and its solutions 

{A-X)(-A* -A) + |B| 2 = (D18) 
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A' = A = - f-7fl ± ^2^2 + |S| 2 - |A| 2 ) (D19) 
irj 7]' V / 



A' = — (-70 ± v/-(0 2 - 1 + 2ntt 2 cos fc) 2 - 2{fl 2 - 1 + 2/d7 2 cos fc)frf} 4 |Q| 2 (cos 2fc - 1)J (D20) 

(We may use the two symbols Q and Q interchangeably in this paper, but it is clear that in a dimensionless equation, 
Q means Q.) 

Appendix E: Dynamic equations for a 2D layered magneto-elastic metamaterial 



15J. It contains two layers (x, y) of 



We consider a two-dimension magneto-elastic metamaterial as in Fig.l of 
split-ring resonators. We will make the following approximations. The mutual magnetic inductances are restricted 
to its counterpart in the other layer, four nearest neighbors in the same layer and the neighbors of the counterpart. 
The mutual mechanic interactions are restricted to its counterpart and its neighbors' counterparts. The mechanic 
interaction coefficients (e.g., the / function in Appendix [2} do not depend on the distance between the split-ring 
resonators. 

The displacement of the SRR at the site (m,n) of the x layer is denoted by x m ^ n . Then 

Kx my n = —fte[(Iy( m ,n)fo + Iy(m+l,n)90 + Iy(m-l,n)90 + Iy{m,n+l)90 + Iy(m,n-\)9o)Ix(m,n)\ ' C^) 

with a similar relation for the y layer. 
Mutual inductances: 



M(x m>n ,x m!n -i) = M(x min ,x mtn+ i) = M{x m , n , x m -i,„) = M{x m , n ,x m+ x, n ) = M) < (E2) 



M(x m . n ,y m . n -i) = M 20 + M 2 i{x m . n - y m ,n-i) (E3) 



M(x m:ni y m ,n+i) = M 20 + M 2 i{x mt n - y m , n +i) (E4) 



M(x m .n, y m -t,n) = M20 + M2\{x m .n - Vm-l,n) 



(E5) 
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M(x m , n , y m +i, n ) = M 2 o + M 2 i(a; m ,„ - y m +i, n ) (E6) 

M{x m , n ,y m<n ) = M10 + M n (i m; „ - y TOi „) (E7) 

The static equation is 

^ =£ + uj' 2 J2 M a Q a - (E8) 

where the summation ^ Q runs over all the neighbors and counterparts of Q x ( m ,n)- 
The dynamic equation under slowly- varying approximation is 

^=£ + w 2 ^M Q Q Q , (E9) 

a 

with a corresponding equation for the y layer. 

The nondimensionalization procedure is similar to that in Appendix [X] with some additional definitions: setting 
/o = 1, Mi,i = 1, and defining g = g /fo, m = M /L, m w = M w /L, m 20 = M 2Q /L, m 2 \ = M 21 /Mu. 
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